Shear zones in granular media: 3D Contact Dynamics simulation 
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Shear zone formation is investigated in slow 3D shear flows. We simulate the linear version of 
the split-bottom shear cell. It is shown that the same type of wide shear zones is achieved in the 
presence as well as in the absence of gravity. We investigate the relaxation of the material towards 
a stationary flow and analyze the stress and the velocity flelds. We provide the functional form 
of the widening of the shear zone inside the bulk. We discuss the growth of the region where the 
material is in critical state. It is found that the growth of the critical zone is responsible for the 
initial transient of the shear zone. 
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I. INTRODUCTION 

Granular materials consist of macroscopic grains which 
have simple interactions and obey the laws of classical 
mechanics. The resulting behavior, however, can be sur- 
prisingly complex. This is shown by a large variety of 
phenomena which, in most cases, present great chal- 
lenges to theoretical descriptions. 

In the present paper we focus on one of the unsettled 
problems of granular media: the quasi-static rheology 
Q. The flow is called quasi-static when inertia effects 
are negligible. This can be achieved by a combination of 
large pressure and low deformation rate. 

An important property observed for this type of flows 
is that stresses become independent of the deformation 
rate. This is remarkable because it is the rate depen- 
dence that is naturally expected to be responsible for the 
rheology. For quasi-static flows the link between stress 
and deformation rate is missing and we lack the consti- 
tutive law which could describe the deformation field in 
the material. 

An experimental setup which is particularly suited to 
provide insight into quasi-static granular flow is the split- 
bottom shear cell. It has recently been the subject of 
many experimental, theoretical and simulational studies 
0, i, i, 0, i, i, [13, im, E d • In experiments the cell 
has cylindrical form called the modified Couette cell Q. 
It is a container whose bottom is divided into a central 
disk and an outer ring. The disk rotates slowly with 
respect to the rest of the container. When sand is filled 
in, it is dragged along by the rotating central bottom disk 
so that a shear zone emerges in the material, where the 
shear deformation is localized. It starts at the perimeter 
of the bottom disk, spreads into the bulk and reaches 
the top surface, if the filling height is not too large. The 
shear zone can be characterized by its central sheet (the 
sheet of the maximum shear rates) and the width of the 
zone around the central sheet. 

Depending on the experimental conditions the behav- 
ior of the central sheet can be quite complicated. Due to 



the cylindrical shearing it gets a nontrivial curved shape 
with decreasing radius towards the top of the system. 
The shape depends strongly on the filling height fS*! . For 
large Hf^w the central sheet even detaches from the top 
and dives entirely into the bulk 0,11, Q forming a cupola- 
like shape. If two materials are used, the central sheet 
can be refracted, when the shear zone leaves one material 
and enters into the other [T2|. AH these effects will be 
avoided hereafter in the paper. We deal with the linear 
version of the split-bottom cell (see later) [l^, [HI , where 
the central sheet remains a vertical plane. Therefore the 
flow becomes simpler and widening of the shear zone can 
be analyzed more easily. 

The width of the shear zone W has been found to be 
an increasing function of the bulk height z and also an 
increasing function of the filling height Hf^w [^, [3, [isj . 
The width at the top (VFtop) grows more slowly than the 
filling height but faster than the square root of -fffin. The 
experimental data suggest that Wtop is approximately a 
power law with exponent 2/3 [1, [1] . 

Although it is a very basic question, where this type of 
rheology comes from, no satisfying description has been 
found so far. There are proposals what the governing 
mechanism of shear zones could be. Some of these the- 
oretical approaches [3, [13, [HI [l3l ^'^^ reassuring in the 
sense that they lead to wide shear zones and, at the same 
time, satisfy the condition of rate independence. One ap- 
proach 0, [12, [Hj is based on the weakest sliding surface 
which fluctuates during the flow, another one [lol . [ll| is 
based on the variation of the effective friction coefficient 
depending on the orientation of the local shear plane. At 
the current stage these models are not very well estab- 
lished and concerning the details they leave many ques- 
tions open. 

It is hard to refine existing models or propose new can- 
didates because not enough details are known about the 
flow. Especially, precise data are needed that are mea- 
sured in the bulk, regarding e.g. velocity and stress flelds. 

With our present study we provide some new details 
about the flow in the split-bottom shear cell. We per- 
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FIG. 1: The linear split-bottom cell. 



form DEM simulations where the velocities and stresses 
are easily accessible in the bulk. Previous experiments 
and simulations were in gravity which leads to an inho- 
mogeneous pressure distribution. We analyze shear zone 
formation also in a zero gravity environment in order to 
clarify the role of gravity. 



II. DESCRIPTION OF THE NUMERICAL 
EXPERIMENT 

In our simulations we examine a linear version of the 
split-bottom shear cell [13, [ll[ shown in Fig. [TJ Here 
the bottom is cut along a straight line. The left and the 
right sides of the boundary move along the y axis in op- 
posite directions both with velocity Ushoar- In y-direction 
periodic boundary conditions are applied. Small grains 
are glued to the side walls and to the bottom in order to 
make their surface rough. 

We have a frictionless piston on the top of the system. 
Its position in z-direction gives the filling height of the 
system. The piston has a compressing force -Fpist on it 
acting in negative z-direction. We use two ways to put 
the system under pressure. Either we apply a large force 
i^pist and set gravity to zero or we use gravity instead 
and put only a weak force on the piston. The role of 
the piston in the former case is to provide the confining 
pressure on the system. In the latter case it keeps only 
the top surface flat and ensures a constant flUing height 
for the whole system. Then the piston has negligible 
effect on the pressure distribution in the bulk which is 
generated essentially by gravity. 

Our simulations are discrete element simulations based 
on the method of contact dynamics [3, . The grains 
are noncohesive, rigid and spherical interacting via fric- 
tional contact forces. The value of the friction coefficient 
is set to 0.2. Throughout this paper every length is mea- 
sured in units of the maximum grain radius. Radii are 
uniformly distributed between 0.8 and 1.0. 

We tested various system sizes. The number of the 
grains N contained by the shear cell varies between 1 000 



and 100 000. The width L^, the length Ly and the filling 
height Hfiw of the systems range from 20 to 240, from 
12.5 to 75, and from 8 to 70, respectively. 

Our simulation corresponds to an experimental situa- 
tion where the grains have density 2400 kg /m? and max- 
imum radius 1 mm. The value of f shear is set to 0.7 
cm/s (unless stated otherwise). The force -Fpist is cho- 
sen proportional to the surface of the piston in order to 
maintain the same pressure. This pressure is 500 N/ni? 
when gravity is switched off. Together with gravity the 
pressure on the piston is set to 25 Njrr? . 

The preparation of the system starts from a gas state 
where grains have random positions. First we compact- 
ify the material with the piston then gravity is switched 
on if needed and the shearing starts. Before measuring 
velocities and stresses we let the system relax in order to 
reach stationary flow. 



III. RESULTS 
A. Orthogonal velocities 

First we examine whether the shear cell generates any 
convection orthogonal to the shear direction y. The com- 
ponents of the coarse-grained velocities Vx , Vy and are 
functions of the coordinates x and z (the coordinate y is 
averaged out). In the present shear cell Vx and Vz would 
vanish for a laminar flow of a Newtonian fluid, however, 
this does not hold a priori for quasi-static flow of granu- 
lar media. One could imagine various kinds of stationary 
flows with non- vanishing convection in x-z-plane, e.g., 
where grains, besides moving in ±?/-direction, slowly rise 
near to the symmetry plane and descend far away from 
it. 

Fig. [2] shows a typical orthogonal velocity field. This 
simulation contains 100000 grains and is performed with- 
out gravity. Before recording velocity data we sheared 
the system for a long time in order to achieve a steady 
state. During this preshearing the system had total shear 
displacement A = 500 (this is the displacement of the two 
sides of the shear cell with respect to each other). The 
velocities shown in Fig. [21a are obtained by an average 
over a further shear displacement 60. This velocity field 
seems to be merely random fiuctuation where vortices 
arise due to mass conservation. The largest velocities are 
located near the symmetry plane of the cell and their 
magnitude is about 200 times smaller than the shear ve- 
locity of the cell. These velocities decrease further if the 
average is taken over larger shear displacements. This is 
shown in Fig. Ob where the shear displacement is 5 times 
larger compared to Fig. [2] a. Thus we find no convection 
orthogonal to the shear direction. If such convection is 
induced by the shearing then it must be at least 4 orders 
of magnitude smaller than the shear velocity itself. 

In the foUowings we will focus on the motion along the 
shear direction y. 
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FIG. 2; The average velocity field in the cross section of the 
cell orthogonal to the shear direction. The time average is 
taken over a shear displacement 60 for (a) and 300 for (b). 
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FIG. 3: The velocity profiles are taken from the same system 
at three different heights z. The line shows an error function 
fit. 
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B. Velocity profiles 

The shear flows found in our simulations are in agree- 
ment with previous experimental and numerical measure- 
ments in the modified Couette cell [1 H, H, H, UH . Our 
new results concern the slow evolution in the outskirts 
of the shear zone, its widening with increasing distance 
from the bottom slit (Sec. C), and the influence of grav- 
ity (Sec. D). In particular it will be shown that gravity 
has surprisingly little effect on the properties of the shear 
zone. 

In order to examine the velocities Vy we divide our sys- 
tem into different slices at constant heights z. z ranges 
from zero to the filling height U'^w. At each height z 
the velocity Vy goes from — Wshoar to fshcar as x is in- 
creased. This transition is very sharp at the bottom, 
where the boundary condition prescribes a step func- 
tion, and broader towards the top of the system. Fig. [3] 
shows the profiles for a system without gravity in several 
heights. 

The velocity profiles can be well fitted with error func- 
tions . Consequently, the shear rate ■^xy as a function 
of a; is a Gaussian curve. We define the width of the shear 
zone W{z) as the square root of the second moment of 
the (normalized) shear rate at height z thus the width of 
the zone equals to the width of the corresponding Gaus- 



FIG. 4: The difference A between the velocity data and the 
fit by error function is shown. The data are recorded in nine 
different systems. Total number of grains used in the simula- 
tions are between 6 000 and 100 000. Minimum and maximum 
heights are 25.6 and 69. In all cases shown here gravity is set 
to zero. Only the fits of the velocity profiles at z — Hfm are 
evaluated. For each system the coordinate x is normalized by 
the width of the shear zone measured in the top layer. 

sian curve. 

The accuracy of the fit by an error function is assessed 
in Figures [4] and [5] In Fig.|4]the deviations between data 
and fit are plotted for several systems. Deviations are 
random and approximately 2% of Ushear near the center. 
Further away from the center the errors become smaller, 
however, systematic deviations can be seen: They are 
positive on the left and negative on the right hand side. 
From Fig. [5] one can conclude that these systematic devi- 
ations are going to vanish, if one lets the simulations run 
longer. The velocity profile approaches the error func- 
tion shape first in the center, but much more slowly in 
the outskirts. If we let the simulations go on, the tail of 
the velocity profile keeps evolving and is getting closer 
and closer to the Gaussian tail of the error function. We 
will come back to the relaxation process in section IIII Fl 

The above properties of the velocity distribution are 
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FIG. 5: Slow time evolution of the tail of the velocity profile 
is shown for a system of 10 000 grains at the filling height 
Hfai — 25.6. For the three different profiles the measure- 
ments are started at shear displacement 200, 800 and 1800. 
Each profile represents an average over an additional shear 
displacement 200. The two lines are fit curves of the velocity 
data: the upper one is a hyperbolic tangent, the lower one is 
the error function. It can be seen that the velocities follow 
the Gaussian tail of the error function the better the larger 
the shear deformation of the sample is. 



in agreement with experimental data [5] which have been 
achieved in gravity in a modified Couette cell, whereas 
the simulation data presented in Fig.s [H [4] and [5] were 
obtained for zero gravity. As we will discuss below, grav- 
ity has indeed only very little influence on the velocity 
profile. 
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FIG. 6: Data collapse of width - bulk height data is achieved 
when rescaled with the maximum width and the filling height, 
respectively. Systems containing different number of grains 
are plotted here. Filling heights ranges from 9 to 69. a) 
without gravity, b) with gravity. 




FIG. 7: The system and its mirrored counterpart above. 



the data follow is a quarter of a circle: 



C. Widening of the shear zone in the bulk 

Previous experimental and numerical studies revealed 
that the shear zone becomes wider as it goes from the 
bottom z = towards the top z = Hf^w. There are some 
experimental and theoretical indications reported in [lol | 
that W[z) may be a power law with exponent between 
0.2 and 0.5. However, we are not aware of any conclu- 
sive experimental data concerning the exact shape of the 
function W{z). It is a crucial question what the func- 
tional form is because it provides a very strong test for 
theories. Such tests are clearly needed as the problem 
of quasi-static flow and shear zone formation is far from 
understood. This field is in the stage of searching for can- 
didates of models in order to gain a better understanding 
of basic phenomena. 

It is a nice feature of computer simulations that one can 
easily access the velocity data also inside the bulk. Based 
on these data we are able to deduce the functional form 
of W{z). We tested many systems with different filling 
heights. It turns out that all the width data collapse if 
plotted in the frame {W{z)/Wtop i z/Hf^\\). The collapse 
of the data can be seen in Fig. [S] The master curve that 



W(z) = Wt 
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Thus we find that the widening of the shear zone starts 
with an exponent 1/2 for small values of z but soon de- 
parts from the power law. W{z) hits the top of the sys- 
tem at a right angle. 

This latter condition of the right angle seems to be 
quite reasonable at least for the case when gravity is 
switched off. The frictionless piston we apply at the top 
exerts no drag force on the material, but only applies 
normal pressure on the system. An equivalent situation 
can be achieved if we take the original system together 
with its mirror image (see Fig. [7]) and at the same time 
we leave the piston away. Then we have a split boundary 
both at the bottom and at the top. The total height of 
the system is then two times the original filling height. 
For symmetry reasons there is no drag force between the 
upper and lower parts of the system which explains the 
equivalence. And again for symmetry reasons the curve 
W{z) must be perpendicular to the plane of the removed 
piston. 

Interestingly the presence or absence of gravity has no 
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FIG. 8: The stress Gzz as the function of depth. Full circles 
are recorded in gravity while open circles correspond to zero 
gravity where the entire pressure is provided by a piston. For 
each case we plot 4 different curves which represent differ- 
ent number of grains filled in the shear cell (A'^: 5000, 10000, 
15000, 20000 and 25000). The fiUing height Jffin changes pro- 
portionally to A^. In gravity ctzz corresponds to hydrostatic 
pressure, ppist is the pressure on the piston in zero gravity. 

influence on the data collapse: the master curve given by 
Eq. [T]is valid for both cases (Fig. [6]). 

D. Role of gravity 

Significant efforts have been made recently to under- 
stand the behavior of wide shear zones. However, all ex- 
perimental, theoretical and numerical studies subjected 
to split bottom shear cells (either linear or cylindrical 
cells) 0, H i, 0, i, i, [13, [O, III investigated shear zones 
under gravity. 

Gravity leads to an inhomogeneous stress distribution 
in the system. Stresses even go to zero as the free sur- 
face of the sample is approached. It is not unplausible 
to imagine that gravity might be responsible for certain 
features of the shear zones (e.g. their widening towards 
the free surface). 

In Fig. [5] we show the stress component a^z at the 
symmetry plane of the shear cell. It can be seen that 
(Jzz is proportional to the depth in case of gravity and 
approximately constant without gravity. 

Recently, Depken et al. [l^, [Hi argued that in the 
quasi-static regime one cannot achieve wide shear zones, 
if the effective friction coefficient /ioff is assumed to be 
constant. For constant /icg the shear zone should localize 
to a thin layer. In their model the widening of the shear 
zone is attributed to the dependence of /icff on the angle 
O between the direction of gravity and the local tangent 
plane of the constant velocity surfaces. (In other words 
they assumed that the frictional properties of the mate- 
rial depend on the orientation of gravity with respect to 
the local sliding plane.) 

Our simulation data do not support the above pic- 
ture. The effective friction /lofr might vary throughout 
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FIG. 9: Width of the shear zone at the top of the systems 
plotted as function of the filling height. Gravity reduces the 
width of the shear zones. Lines are not fits, they only show 
the slope of exponent 2/3. 



the shear zone, however, the direction of gravity does not 
seem to play any important role here. As we discussed 
in previous sections the velocity profiles are qualitatively 
the same no matter whether gravity is present or not. 
Surprisingly, body forces or decreasing pressure are not 
needed for wide shear zones. 

In fact, shear zones exhibit even larger width when 
gravity is switched off. The width gets larger by a factor 
1.2 ±0.1. This can bee seen in Fig. [5] where the top width 
of the shear zone for various filling heights is plotted. 

Why gravity contracts the shear zone can be explained 
by the pressure distribution. At the bottom the pressure, 
and thus also the frictional forces, are much stronger than 
at the top of the system. Therefore the rheology is dom- 
inated by the lower part of the system. Compared to 
the case of homogeneous pressure this acts as if the sys- 
tem had effectively smaller filling height which leads to 
a smaller width of the shear zone. The contracting effect 
of gravity is also discussed in fl6] based on the principle 
of minimum energy dissipation. 

The experimental data of Fenistein et al. [1, Q showed 
that Wtop(^^fiii) is approximately a power law with expo- 
nent 2/3. W^top(^^fiii) found in our simulations is shown 
in a log- log plot in Fig. [51 where the exponent 2/3 is also 
indicated for comparison. The data follow approximately 
the experimental behavior. For a precision value of the 
exponent or a discussion of deviations from a power law 
better statistics is needed, however. 



E. Influence of additional parameters 

It was assumed so far that the flow is quasi-static, i.e. 
the shear velocity tighear is small enough that no rate- 
dependence is observed in the behavior of the shear zone. 
We also intended to choose the width L^, and length Ly 
of the systems large enough in order to exclude their in- 
fluence on the flow. Furthermore, it was assumed that 
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FIG. 10: Robustness of the width of the shear zone, a) The 
effect of the system width L^c on the top width of the shear 
zone, b) The width of the shear zone is plotted as function of 
the bulk height z for three different lengths Ly of the shear 
cell, c) Relaxation of the top width as function of the shear 
displacement A. d) The influence of the shear velocity Hahcar 
on the curve W{z). FuU circles denote five times larger shear 
velocity than the open circles. 



we start collecting data. Fig. [TOlc shows that the tran- 
sient period is indeed completed during the preshearing: 
The width of the shear zone becomes independent of the 
displacement A. The origin of the initial transient is dis- 
cussed in Sec. IIIIFI 

Due to the larger time resolution applied here veloc- 
ity fluctuations become also larger. In order to improve 
the statistics 20 parallel simulations were taken. In each 
case the system had the same macroscopic parameters 
as our reference system. In Fig. llOl c the time evolution 
of Wtop was achieved by ensemble average over the 20 
simulations. 

Next we turn to Wshcar- If pressure conditions are the 
same it depends on the speed of driving whether the 
flow is quasi-static. For relative large driving speeds, 
but still in the dense flow regime, inertia effects come 
into play which weaken shear localization and results in 
wider shear zones. Then the width of the shear zone can 
be made smaller by lowering the driving speed. But if 
the driving gets slow enough the rheology becomes inde- 
pendent of the driving rate and the shear zone reaches 
its minimum width. We test the effect of the driving 
rate in our reference system by taking five times smaller 
shear velocity. The new bulk-width of the shear zone is 
compared with the original one in Fig. [TOld. The data 
show no further decrease of the width. The width for 
the reduced driving is even slightly larger due to ran- 
dom fluctuations. Within the accuracy of our numerical 
measurement the two curves can be regarded as equal. 



the analysis of the shear zone was taken after initial tran- 
sients in the stationary flow regime, i.e. the total shear 
displacement A was large enough to ensure complete re- 
laxation. 

In this section we check the influence of the parameters 

^ and t'shcar 

to show that they are chosen prop- 
erly and do not alter the properties of the shear zone. For 
this purpose we take one of the previous samples as refer- 
ence system and demonstrate the role of the four param- 
eters there. The reference system contains N = 10 000 
grains, has the size = 80, Ly = 25, Hf^\\ = 25.6 (this 
system appeared already in Figures HI El [6] and (H) . 

First we vary the width of the system and the num- 
ber of grains TV proportionally to in order to maintain 
the same filling height. Other parameters are constant. 
The effect on the width of the shear zone can bee seen in 
Fig. [Tola. Around the reference point = 80 the width 
remains constant. The shear zone "feels" the effect of the 
side walls only if drops below 40 where the shear zone 
can be strongly contracted by lowering L^. 

For the parameter Ly we test values 12.5, 25 and 50 
(again N is changed proportionally) while other param- 
eters are constant. Fig. [TOlb shows that the width of 
the shear zone remains essentially the same (changes are 
small and without systematics) . 

In order to avoid initial transients all the systems pre- 
sented in this paper undergo a preshearing over a total 
displacement of more than A > 60, and only after that 



F. Shear zone versus critical zone 

The shear zone is the region where the major part of 
the shear deformation takes place. It is described by 
the function W{z). The material, however, is not solid 
outside the shear zone either. Split-bottom cells fluidize 
the material everywhere, however the shear rate becomes 
many orders of magnitude smaller far from the shear 
zone. 

In this section we would like to discuss the concept 
of the critical state. It has not got any attention so far 
in the context of the split-bottom shear cells, although, 
it leads to the emergence of a relevant and new type of 
zone. 

It is known that the mechanical properties of granular 
media are influenced strongly by the preparation. If one 
starts shearing a packing the behavior can be different 
depending on the initial state (density, structure of the 
contact network, etc): It can lead to different stress re- 
sponses, effective frictions, dilation or contraction, etc. 
However, if the material experiences large enough local 
strain it reaches a unique state regardless of the prepara- 
tion history. This is the critical state [13, H [M 0, HH 
where the material organizes and maintains its micro- 
scopic inner structure on shearing. After the critical state 
has been reached unlimited shear deformation can occur 
without changes of stresses or density. The characteris- 
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tic deformation scale needed to erase the memory of the 
material and reach the critical state is typically around 
7 — 0.2, where 7 is the cumulative shear strain. 

Before the shearing starts in the split-bottom cell the 
fabric of the material reflects the direction of the ini- 
tial compression or gravity. With the shear deformation 
this structure is destroyed and new contacts are created 
against the direction of the shear. This gives rise to strain 
hardening: the resistance of the material against shear is 
increased. This does not happen simultaneously all over 
the sample. When regions in the middle of the shear cell 
are already in the critical state, regions far away can still 
be frozen in the initial configuration. 

At the beginning, the zone of the critical state starts 
growing from the split line at the bottom. It reaches 
quickly the top of the system and also spreads sideways. 
As the shear rate is very small far from the symmetry 
plane, the growth of the critical zone becomes extremely 
slow here. In that sense a steady state can not be reached 
in the whole system. One expects the flow to become 
stationary only inside the critical zone. 

The growth of the critical zone comes into sight in 
Fig. [5l The velocity of a given point in the system reaches 
its final value only after the region becomes critical. The 
position where the velocity data depart from the sta- 
tionary curve mark the border of the critical zone. The 
different velocity profiles recorded in different stages of 
the simulation show how the width of the critical zone is 
increasing with time. It can be attributed to the strain 
hardening why the shear rate is reduced at a given posi- 
tion after it becomes critical. 

The reduction of the shear rate due to strain harden- 
ing also explains the relaxation of the shear zone that is 
shown in Fig. [TUlc. In the early stage of the simulation 
where the critical zone is smaller than the shear zone the 
shear rate is slightly enhanced outside the critical region. 
This makes the shear zone wider a little bit. After the 
whole shear zone becomes critical this additional widen- 
ing effect ceases, and Wtop is reduced to its final value. 

According to this interpretation the initial transient of 
the shear zone is due the time evolution of the critical 
region. One can estimate the shear displacement A that 
corresponds to the transient period by matching the size 
of the stationary shear zone and the growing critical zone. 
Using the stationary velocity profile and the cutoff shear 
strain 7 = 0.2 one gets A = 7.7 for the transient. This 
is in excellent agreement with the relaxation observed in 
Fig. Mc. 

The growth of the critical zone and the transient of the 
shear zone can be reproduced within the framework of a 
simple lattice model [l^ developed for quasi-static shear 
flows. 

The presence of the critical zone can be observed also 
in the stress field. The next section is subjected to this 
question. 



G. Stresses 

As the average local velocities have only a y-component 
ant the system is translational invariant in j/-direction, 
the local strain rate in the (x, y, z)-frame has the form 
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By means of a rotation around the y-axis a new local 
frame (u, y, v) can be chosen in which the strain rate is 







dVy_ 

du 







^00 




(3) 



If initial conditions are forgotten in the critical state 
(also called steady state flow), stress and fabric tensors 
are expected to have the same principal axes as the strain 
rate tensor. Then the stress tensor must have the form 



m 













T 







^yy 


<7yv 






P 







^vy 






^0 





p 



(4) 



Depken et al. tested stresses in the linear split bot- 
tom cell by soft particle molecular dynamics simulations. 
They found the expected behavior in the middle part of 
the cell where stress and strain tensor were co-linear and 
stresses took the form (|4]). 

However, if the material still remembers its initial 
structure the alignment of stress and strain is not neces- 
sarily valid. Therefore we determine a second local frame 
{u',y,v') from the condition ayyi — a^iy — and eval- 
uate the angle a between the two directions (u, v) and 
{u',v') in the (a;, 2:)-plane. The lower part of Fig. [Tl] 
shows indeed that the angle a approaches for large 
shear deformations 7, while the principal directions of 
stress and strain rate tensor differ during the transient. 
Correspondingly, au'v' approaches for large shear de- 
formations (see Fig. [T2|) , as predicted by EqHl 

Our numerical test was based on another method than 
the one used by Depken et al. [HI, the contact dynamics 
algorithm [3, [3], and we used slightly different con- 
ditions (zero gravity, piston). Nonetheless, we found the 
same behavior for regions where the material experienced 
large shear deformation. Here stress and strain tensors 
align and the stress corresponds to the reduced form in 
Eq.H 

Stress data recorded in a system of 100 000 grains with 
total shear displacement 820 are presented here as the 
function of the cumulative local shear strain 7. Fig. [TT] 
shows stress ratios cjyy/au'u' and a^iyi /<Ju'u' ■ In the crit- 
ical zone CTu'm' and indeed the same and the value 
ayiyi is about 10% smaller. 
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FIG. 11: The effect of the local shear strain 7. Crosses and 
circles show ratios of normal stress components Oyy/u^^^i and 
o'v'v' / <^u'u' , respectively. The parameter a (triangles) indi- 
cate the angle between the local shear stress and shear strain. 
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FIG. 12: The remaining shear stress a^i^i vanishes for large 
shear deformations. The open circles and dots correspond to 
measurements started at A = 285 and 700, respectively. In 
both cases the measurement lasted over a period A = 120. 

IV. CONCLUSIONS 

We studied shear flows in a linear split-bottom cell 
by means of computer simulations. The formation of 



wide shear zones was analyzed in the presence and in the 
absence of gravity. In the former case pressure scales with 
depth, in the latter case it is approximately constant. 
However, in both cases the same type of wide shear zones 
emerge. 

We showed that the widening of the shear zone in the 
bulk can be described by one master curve which holds 
for various sizes and pressure conditions. The shape of 
the widening function is a quarter of a circle and not a 
power law as was suggested before. We hope that this 
result will promote the development of the proper con- 
tinuum theory for quasi-static flows. 

We analyzed the persistent growth of the critical zone 
and its effect on the rheology. It influences the velocity 
field and the stresses, and it causes a transient of the 
shear zone at the beginning of the shear test. 

It was shown that the form of the stress tensor becomes 
simpler with increasing shear strain 7. The region, where 
the stress tensor takes the reduced form (|4]), was also an- 
alyzed by Depken et al. [ll| . The authors found that this 
region can be best characterized by the inertial parame- 
ter / [s'l which is defined to be proportional to the shear 
ratej7 and to the inverse pressure. As it was pointed out 
in 11] it is not clear how the emergence of the inertial 
number can be reconciled with the rate independence of 
quasi-static flows. It is the task of future work to clarify 
the question what influence the parameters 7 and / have 
on the stress tensor in case of slow deformations. 
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